//这里在全部染完之后进行重染色，从最大的颜色进行重染色
//这个程序执行理论最优的图染色算法
#include<stdio.h>
#include<stdlib.h>
#include<string.h>
#include<cusparse.h>
#include<sys/time.h>

#define A_row_size 262144

__global__ void init_need_to_color_device(int* need_to_color_device, int size){
    //每个线程为一个位置附上一个数字
    int thread_id = threadIdx.x + blockIdx.x * blockDim.x;

    //附上一个数字
    if(thread_id < size){
        need_to_color_device[thread_id] = thread_id;
    }
}

//用一个函数来得到染色的最大节点
//如果对应颜色编号有颜色，那么我们就将color_arr的对应为置为1
__global__ void get_color_count(int* coloring, int row_num, int* color_arr){
    //计算当前节点的线程号
    int thread_id = threadIdx.x + blockIdx.x * blockDim.x;

    int color = coloring[thread_id];

    //如果颜色没有登记，那么我们就登记颜色
    if(color_arr[color] == 0){
        color_arr[color] = 1;
    }
}

//对某一个颜色进行重染色
__global__ void Re_color(int* A_device, size_t pitch,int* element_num_each_row_device, 
int row_num, int* coloring, int target_color){
    //计算当前线程号，每个线程负责一个节点
    int thread_id = threadIdx.x + blockIdx.x * blockDim.x;

    if(thread_id < row_num && coloring[thread_id] == target_color){
        //这里进行重染色
        int target_point = thread_id;

        //当前节点的邻居节点的数量
        int neighbor_num = element_num_each_row_device[target_point];

        //首先查看所有邻居节点的染色情况
        
        //使用一个64位数字来保存周围节点已经染的颜色
        //在某一位上是1代表这一位所代表的颜色已经染色了
        //一开始初始化为0，那么说明已有的颜色占用初始化为0
        long neighbor_colors_bit = 0;

        //现在遍历所有的节点
        int i;

        //邻居节点颜色的临时变量
        int neighbor_color;

        //邻居节点的编号
        int neighbor;

        //将矩阵A中的某一行拿出来
        int* A_device_row = (int*)((char*)A_device + target_point * pitch); 

        //存某一个邻居节点的颜色的bit
        long color_bit;

        for(i = 0; i < neighbor_num; i++){
            neighbor = A_device_row[i];

            //通过邻居节点的编号获取邻居节点的颜色，没有染色那么邻居节点就是0
            neighbor_color = coloring[neighbor];

            //如果邻居节点已经染色，那么我们就更新当前节点不能染的颜色
            if(neighbor_color != 0 && neighbor != target_point){
                //已经染色，左移
                color_bit = 1 << (neighbor_color - 1);

                //执行或运算，记录当前那些颜色已经被占用
                neighbor_colors_bit = neighbor_colors_bit | color_bit;
            }
        }

        //这里已经记录了周围节点的所有颜色，然后再进行移位计算，找出可以染色的最小颜色编号
        //我们可以通过mod2来获取最低位的bit位，我们一直移位，直到最低位的bit位为0
        
        //下面这个数字存储了当前节点需要染的颜色，一开始默认是1
        //随着移位操作的不断进行，这个数字会越来越大
        int target_point_color = 1;

        while(neighbor_colors_bit % 2 == 1){
            //当前颜色已经被染色了，移位
            neighbor_colors_bit = neighbor_colors_bit >> 1;
            
            //下一个要查看的颜色
            target_point_color++;

            if(target_point_color > 63){
                printf("出错了，不用那么多颜色\n");
            }
        }

        //当前节点的颜色已经拿到了
        coloring[target_point] = target_point_color;
    }
}

//不进行bit压缩的方案
__global__ void AssignColors_no_compress(int* A_device, size_t pitch,int* element_num_each_row_device, 
int row_num, int* need_color_num, int* coloring, int *need_to_color){
    //计算当前线程号，每个线程负责一个节点
    int thread_id = threadIdx.x + blockIdx.x * blockDim.x;

    //排除掉不存在的点
    if(thread_id < *need_color_num){
        //获取当前节点的编号
        int target_point = need_to_color[thread_id];

        //当前节点的邻居节点的数量
        int neighbor_num = element_num_each_row_device[target_point];

        //首先查看所有邻居节点的染色情况
        
        //使用一个64位数字来保存周围节点已经染的颜色
        //在某一位上是1代表这一位所代表的颜色已经染色了
        //一开始初始化为0，那么说明已有的颜色占用初始化为0
        long neighbor_colors_bit = 0;

        //申请一个数组来存储周围节点的情况
        int can_not[27];

        //现在遍历所有的节点
        int i;

        for(i = 0; i < 27; i++){
            can_not[i] = 0;
        }

        //邻居节点颜色的临时变量
        int neighbor_color;

        //邻居节点的编号
        int neighbor;

        //将矩阵A中的某一行拿出来
        int* A_device_row = (int*)((char*)A_device + target_point * pitch); 

        //存某一个邻居节点的颜色的bit
        long color_bit;

        for(i = 0; i < neighbor_num; i++){
            neighbor = A_device_row[i];

            //通过邻居节点的编号获取邻居节点的颜色，没有染色那么邻居节点就是0
            neighbor_color = coloring[neighbor];

            //如果邻居节点已经染色，那么我们就更新当前节点不能染的颜色
            if(neighbor_color != 0){
                //已经染色
                can_not[neighbor_color] = 1;
            }
        }

        //这里已经记录了周围节点的所有颜色，然后再进行移位计算，找出可以染色的最小颜色编号
        //我们可以通过mod2来获取最低位的bit位，我们一直移位，直到最低位的bit位为0
        
        //下面这个数字存储了当前节点需要染的颜色，一开始默认是1
        //随着移位操作的不断进行，这个数字会越来越大
        int target_point_color = 1;

        while(can_not[target_point_color] == 1){
            target_point_color++;

            if(target_point_color > 63){
                printf("出错了，不用那么多颜色\n");
            }
        }

        //当前节点的颜色已经拿到了
        coloring[target_point] = target_point_color;

        //对于已染色节点的数量，在检测冲突的时候更新
    }
}


//根据周围节点染色的颜色来染一个序号最小的颜色
//传入参数是没有染色的节点列表，经过CSR压缩的图邻接矩阵A
//每个节点的染色情况数组
//我们要使用一个long数字来存储这个线程所负责的节点的邻居节点的颜色
//int* A_device, size_t pitch,int* element_num_each_row_device, int row_num这几个形参代表了图的邻接矩阵A
//int* need_color_num代表了还没有染色的节点的数量
//int* coloring代表了每个点染色的情况
//然后我们需要一个数组来存储没有染色的点的列表
__global__ void AssignColors(int* A_device, size_t pitch,int* element_num_each_row_device, 
int row_num, int* need_color_num, int* coloring_old, int* coloring,int *need_to_color){
        
    //计算当前线程号，每个线程负责一个节点
    int thread_id = threadIdx.x + blockIdx.x * blockDim.x;
    
    //排除掉不存在的点
    if(thread_id < *need_color_num){
        //获取当前节点的编号
        int target_point = need_to_color[thread_id];

        //当前节点的邻居节点的数量
        int neighbor_num = element_num_each_row_device[target_point];

        //首先查看所有邻居节点的染色情况
        
        //使用一个64位数字来保存周围节点已经染的颜色
        //在某一位上是1代表这一位所代表的颜色已经染色了
        //一开始初始化为0，那么说明已有的颜色占用初始化为0
        long neighbor_colors_bit = 0;

        //现在遍历所有的节点
        int i;

        //邻居节点颜色的临时变量
        int neighbor_color;

        //邻居节点的编号
        int neighbor;

        //将矩阵A中的某一行拿出来
        int* A_device_row = (int*)((char*)A_device + target_point * pitch); 

        //存某一个邻居节点的颜色的bit
        long color_bit;

        for(i = 0; i < neighbor_num; i++){
            neighbor = A_device_row[i];

            //通过邻居节点的编号获取邻居节点的颜色，没有染色那么邻居节点就是0
            neighbor_color = coloring_old[neighbor];

            //如果邻居节点已经染色，那么我们就更新当前节点不能染的颜色
            if(neighbor_color != 0){
                //已经染色，左移
                color_bit = 1 << (neighbor_color - 1);

                //执行或运算，记录当前那些颜色已经被占用
                neighbor_colors_bit = neighbor_colors_bit | color_bit;
            }
        }

        //这里已经记录了周围节点的所有颜色，然后再进行移位计算，找出可以染色的最小颜色编号
        //我们可以通过mod2来获取最低位的bit位，我们一直移位，直到最低位的bit位为0
        
        //下面这个数字存储了当前节点需要染的颜色，一开始默认是1
        //随着移位操作的不断进行，这个数字会越来越大
        int target_point_color = 1;

        while(neighbor_colors_bit % 2 == 1){
            //当前颜色已经被染色了，移位
            neighbor_colors_bit = neighbor_colors_bit >> 1;
            
            //下一个要查看的颜色
            target_point_color++;

            if(target_point_color > 63){
                printf("出错了，不用那么多颜色\n");
            }
        }

        //当前节点的颜色已经拿到了
        coloring[target_point] = target_point_color;

        //对于已染色节点的数量，在检测冲突的时候更新
    }
}

//这里执行纠错工作
//传入的参数几乎没有区别，主要就是更新need_color_num和need_to_color两个变量，已经将一部分coloring的点取消掉
__global__ void JudgeConflict(int* A_device, size_t pitch,int* element_num_each_row_device, 
int row_num, int* need_color_num, int* coloring, int *need_to_color, int *need_to_color_new,
int* need_color_num_new){
    //用一个线程扫描所有的节点，然后更新need_to_color
    //计算当前线程号，每个线程负责上次受到染色的节点，然后得出一个新的需要染色节点的列表
    //我们只需要计算当前节点是不是需要进行下一轮染色就可以了
    int thread_id = threadIdx.x + blockIdx.x * blockDim.x;

    if(thread_id < *need_color_num){
        //获取当前节点的编号
        int target_point = need_to_color[thread_id];

        //当前节点的邻居节点的数量
        int neighbor_num = element_num_each_row_device[target_point];

        //获取邻居节点的编号
        int neighbor;

        //获取邻居节点的颜色的临时变量
        int neighbor_color;

        //获取当前节点的颜色
        int target_point_color = coloring[target_point];

        //将矩阵A中的某一行拿出来
        int* A_device_row = (int*)((char*)A_device + target_point * pitch); 

        //遍历周围所有节点的颜色
        int i;
        for(i = 0; i < neighbor_num; i++){
            //获取节点编号
            neighbor = A_device_row[i];

            //获取节点颜色
            neighbor_color = coloring[neighbor];

            //如果当前节点与邻接节点颜色一样，并且需要比邻接节点大，那么我们就要将当前节点加入下一轮染色的序列
            if(neighbor_color == target_point_color && target_point > neighbor){
                //这里说明当前节点要放弃颜色了，需要进行原子性操作来增加need_color_num_new
                int insert_index = atomicAdd(need_color_num_new, 1);

                //插入新的元素
                need_to_color_new[insert_index] = target_point;

                //删除自己的颜色
                coloring[target_point] = 0;

                //退出循环
                break;
            }
        }
    }
}

//检查颜色是不是打印正确的内核函数
//传入矩阵A以及每个节点的染色情况，看看某一个节点的邻居节点的颜色是不是冲突了
__global__ void check_graph_coloring(int* A_device, size_t pitch, int* element_num_each_row_device, 
    int row_num, int* coloring){
    //遍历每一行
    int thread_id = threadIdx.x + blockIdx.x * blockDim.x;
    
    if(thread_id < row_num){
        //将一行取出
        int* A_device_row = (int*)((char*)A_device + thread_id * pitch);

        int col_num = element_num_each_row_device[thread_id];

        int neighbor_col;

        int i;

        //获取当前节点的颜色
        int color = coloring[thread_id];

        if(color == 0){
            printf("出现了错误，有节点没有染色\n");
        }

        for(i = 0; i < col_num; i++){
            //查看当前节点的列号
            neighbor_col = A_device_row[i];

            if(coloring[neighbor_col] != 0 && color != 0 && color == coloring[neighbor_col] && neighbor_col != thread_id){
                //这里说明染色冲突
                printf("出现错了，%d与%d冲突，冲突颜色%d\n", thread_id, neighbor_col, color);
            }
        }
    }
}


int main(){
    //我们创造一个二维数组，这个二维数组一共4096行，27列，每一行就代表了图的节点数量
    //每一列就代表了就是一个点，一行中的不同元素代表了这个节点和哪些行有连接。
    //创造4096个指针
    //把A改成堆区
    int** A = 0;

    cudaMallocHost((void**)&A, A_row_size * sizeof(int *));

    //为每个指针创造27个空间
    //遍历变量
    printf("申请空间\n");

    int i;
    for(i = 0; i < A_row_size; i++){
        cudaMallocHost((void **)&A[i], 27 * sizeof(int));
    }

    //用一个堆区数组来记录每一行元素的数量
    int* element_num_each_row = 0;

    cudaMallocHost((void **)&element_num_each_row, A_row_size * sizeof(int));

    cudaMemset(element_num_each_row, 0, A_row_size * sizeof(int));
    
    //我们将数据写入。
    
    //从文件中读入数据
    FILE* Graph_File;

    //打开文件
    Graph_File = fopen("../data/A_Matrix_col_1.txt", "r");

    //用一个缓冲区来获取读取的内容
    char read_line_buf[1024];

    //用一个变量来遍历A的每一行
    int A_row_index = 0;

    printf("读取文件\n");

    fgets(read_line_buf, 1024, Graph_File);

    while (!feof(Graph_File)) {
        
        //用一个变量来定位A中每一行的每一个元素
        int A_col_index = 0;
        
        //分割字符串，这里是每一个子串
        char* part_char;

        part_char = strtok(read_line_buf, " ");


        while(part_char != NULL && A_col_index < 27){

            //第一个元素
            A[A_row_index][A_col_index] = atoi(part_char);

            // printf("%d ", A[A_row_index][A_col_index]);
            
            A_col_index++;
            element_num_each_row[A_row_index]++;

            part_char = strtok(NULL, " ");
        }

        //先读，才能查到eof
        fgets(read_line_buf, 1024, Graph_File);

        A_row_index++;
    }

    fclose(Graph_File);

    //首先先使用一个内核函数来为每个节点染色
    //这个函数执行两个事情，一个是查看周围节点的染色情况，一个是根据周围节点的染色情况来决定自己应该染的颜色
    //设备端A矩阵指针
    int* A_device = 0;
    //每一行对齐之后的大小
    size_t pitch;
    
    printf("申请设备端的A指针\n");
    //申请设备端的A指针
    cudaMallocPitch((void**)&A_device, &pitch, 27 * sizeof(int), A_row_size);

    //将数组A的数据拷贝给设备端，值的一提的是，数组A的每一行在存储上不是连续的，所以我们使用memcpy2D来处理
    //我们所能用的方式就是将A矩阵一行一行地传给设备端
    printf("pitch = %d\n", pitch);

    for(i = 0; i < A_row_size; i++){
        //因为pitch是字节数量，如果想让A_device这种int指针做地址偏移的话，需要将A_device的指针类型变化一下。
        cudaMemcpy((char *)A_device + i * pitch, A[i], 
            element_num_each_row[i] * sizeof(int), cudaMemcpyHostToDevice);
    }
    cudaThreadSynchronize();


    printf("将数组A的数据拷贝给设备端完毕\n");

    //然后将每一行多少个元素的数组拷贝到GPU端
    int* element_num_each_row_device;

    cudaMalloc((void **)&element_num_each_row_device, A_row_size * sizeof(int));

    //将数据拷贝到设备端
    cudaMemcpy((void *)element_num_each_row_device, element_num_each_row, A_row_size * sizeof(int), cudaMemcpyHostToDevice);

    printf("然后将每一行多少个元素的数组拷贝到GPU端\n");

    //创建一个变量来保存还没有染色的节点的数量
    int need_color_num = A_row_size;

    int* need_color_num_device;
    int* need_color_num_new_device;

    //将这个变量传到显存上
    cudaMalloc((void**)&need_color_num_device, sizeof(int));
    cudaMalloc((void**)&need_color_num_new_device, sizeof(int));

    cudaMemcpy(need_color_num_device, &need_color_num, sizeof(int), cudaMemcpyHostToDevice);

    //创建一个数组来表达每个节点所染的颜色
    int* coloring = 0;
    
    cudaMallocHost((void **)&coloring, sizeof(int) * A_row_size);

    int* coloring_device = 0;
    //在创建一个一样的数组来存储上一轮染色的结果
    int* coloring_device_old = 0;

    cudaMalloc((void **)&coloring_device, sizeof(int) * A_row_size);
    cudaMalloc((void **)&coloring_device_old, sizeof(int) * A_row_size);
    //将整块显存初始化为0
    cudaMemset(coloring_device, 0, sizeof(int) * A_row_size);
    cudaMemset(coloring_device_old, 0, sizeof(int) * A_row_size);

    //创建颜色数量
    int color_num = 0;

    //创建线程结构，线性网格，每个block64个线程，每次迭代之后都需要更新
    dim3 block_dimension(384);

    dim3 grid_dimension(A_row_size / 384 + 1);

    //这里存储还没有染色的节点
    //我们只在GPU端申请空间
    int* need_to_color;
    int* need_to_color_device;
    int* need_to_color_new_device;

    //申请空间，空间大小和A和行数一样
    //这个数组的实际大小就是总行数-已经染色的节点数量
    cudaMallocHost((void**)&need_to_color, A_row_size * sizeof(int));
    cudaMalloc((void**)&need_to_color_device, A_row_size * sizeof(int));
    cudaMalloc((void**)&need_to_color_new_device, A_row_size * sizeof(int));

    //我们使用一个gpu函数来初始化这个数组
    init_need_to_color_device<<<grid_dimension, block_dimension>>>(need_to_color_device, A_row_size);

    // //将数组传到cpu端，然后打印
    // int* need_to_color;

    // cudaMallocHost((void **)&need_to_color, A_row_size * sizeof(int));

    // //将数据拷贝到cpu端
    // cudaMemcpy(need_to_color, need_to_color_device, A_row_size * sizeof(int), cudaMemcpyDeviceToHost);

    // for(i = 0; i < A_row_size; i++){
    //     printf("%d ", need_to_color[i]);
    // }

    // printf("\n ", need_to_color[i]);
    
    //开始执行染色函数
    //第一次执行应该染成颜色1
    //AssignColors(int* A_device, size_t pitch,int* element_num_each_row_device, 
    //int row_num, int* count_device, int* coloring, int *need_to_color)
    //计时开始
    struct timeval starttime,endtime;
    
    double timeuse;

    gettimeofday(&starttime,NULL);

    while(need_color_num != 0){
        printf("染色，need_color_num：%d\n", need_color_num);

        cudaMemcpy(coloring_device_old, coloring_device, A_row_size * sizeof(int), cudaMemcpyDeviceToDevice);

        AssignColors<<<need_color_num / 384 + 1, block_dimension>>>(A_device, pitch, element_num_each_row_device, A_row_size, 
            need_color_num_device, coloring_device_old, coloring_device, need_to_color_device);
        
        //然后进行纠错
        //void JudgeConflict(int* A_device, size_t pitch,int* element_num_each_row_device, 
        //int row_num, int* need_color_num, int* coloring, int *need_to_color, int *need_to_color_new,
        // int* need_color_num_new)
        //初始化need_color_num_new_device
        cudaMemset(need_color_num_new_device, 0, sizeof(int));

        JudgeConflict<<<need_color_num / 384 + 1, block_dimension>>>(A_device, pitch, element_num_each_row_device, A_row_size,
            need_color_num_device, coloring_device, need_to_color_device, need_to_color_new_device, need_color_num_new_device);
        
        //更新新的待处理节点和待处理节点数量
        cudaMemcpy(need_color_num_device, need_color_num_new_device, sizeof(int), cudaMemcpyDeviceToDevice);
        cudaMemcpy(need_to_color_device, need_to_color_new_device, sizeof(int) * A_row_size, cudaMemcpyDeviceToDevice);
        
        //将need_color_num拷贝回来
        cudaMemcpy(&need_color_num, need_color_num_device, sizeof(int), cudaMemcpyDeviceToHost);

        //将数据拷贝到cpu端
        // cudaMemcpy(need_to_color, need_to_color_device, A_row_size * sizeof(int), cudaMemcpyDeviceToHost);

        // for(i = 0; i < need_color_num; i++){
        //     printf("%d ", need_to_color[i]);
        // }

        // printf("\n\n ", need_to_color[i]);

        //更新线程结构
    }

    

    //将染色情况传回，查看染色情况
    cudaMemcpy(coloring, coloring_device, A_row_size * sizeof(int), cudaMemcpyDeviceToHost);

    // 打印染色情况
    // for(i = 0; i < A_row_size; i++){
    //     printf("%d ", coloring[i]);
    // }

    // printf("\n");

    //展开自检
    //check_graph_coloring(int* A_device, size_t pitch, int* element_num_each_row_device, int row_num, int* coloring)
    check_graph_coloring<<<A_row_size / 384 + 1, block_dimension>>>(A_device, pitch, element_num_each_row_device, A_row_size, coloring_device);
    
    //申请一个数组
    int* color_arr;
    int* color_arr_device;

    cudaMalloc((void **)&color_arr_device, 50 * sizeof(int));
    cudaMallocHost((void **)&color_arr, 50 * sizeof(int));

    //初始化color_arr_device
    cudaMemset(color_arr_device, 0, 50 * sizeof(int));
    
    //get_color_count(int* coloring, int row_num, int* color_arr)
    get_color_count<<<A_row_size / 384 + 1, block_dimension>>>(coloring_device, A_row_size, color_arr_device);

    //将结果传回来
    cudaMemcpy(color_arr, color_arr_device, 50 * sizeof(int), cudaMemcpyDeviceToHost);

    color_num = 0;

    for(i = 0 ; i < 50; i++){
        printf("%d ", color_arr[i]);
        if(color_arr[i] == 1){
            color_num++;
        }
    }

    printf("\n颜色的数量：%d\n", color_num);

    //这里执行重染色
    for(i = 2; i <= color_num ; i++){
        //void Re_color(int* A_device, size_t pitch,int* element_num_each_row_device, 
        //int row_num, int* coloring, int target_color)
        Re_color<<<A_row_size / 384 + 1, block_dimension>>>(A_device, pitch, element_num_each_row_device,
            A_row_size, coloring_device, i);
    }

    //计时结束
    gettimeofday(&endtime,NULL);

    timeuse=(1000000*(endtime.tv_sec-starttime.tv_sec)+endtime.tv_usec-starttime.tv_usec)/1000;
    
    printf("使用的时间:%fms\n", timeuse);

    //再次查看颜色数量
    //初始化color_arr_device
    cudaMemset(color_arr_device, 0, 50 * sizeof(int));

    get_color_count<<<A_row_size / 384 + 1, block_dimension>>>(coloring_device, A_row_size, color_arr_device);

    //将结果传回来
    cudaMemcpy(color_arr, color_arr_device, 50 * sizeof(int), cudaMemcpyDeviceToHost);

    color_num = 0;

    for(i = 0 ; i < 50; i++){
        printf("%d ", color_arr[i]);
        if(color_arr[i] == 1){
            color_num++;
        }
    }

    printf("\n颜色的数量：%d\n", color_num);

    //析构need_to_color_device
    cudaFree(need_to_color_device);

    //析构need_to_color
    cudaFreeHost(need_to_color);

    //析构coloring_device
    cudaFree(coloring_device);

    for(i = 0; i < A_row_size; i++){
        cudaFreeHost(A[i]);
    }

    //析构矩阵A
    cudaFreeHost(A);

    //析构coloring
    cudaFreeHost(coloring);

    //析构coloring_device
    cudaFree(coloring_device);
    cudaFree(coloring_device_old);

    //析构count_device
    cudaFree(need_color_num_device);

    //析构need_color_num_new_device
    cudaFree(need_color_num_new_device);

    //析构element_num_each_row_device
    cudaFree(element_num_each_row_device);

    //析构element_num_each_row
    cudaFreeHost(element_num_each_row);

    //析构need_color_num_new_device
    cudaFree(need_color_num_new_device);

}



